Event-chain Monte Carlo algorithms for hard-sphere systems * 
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In this paper we present the event-chain algorithms, which are fast Markov-chain Monte Carlo 
methods for hard spheres and related systems. In a single move of these rejection- free methods, an 
arbitrarily long chain of particles is displaced, and long-range coherent motion can be induced. Nu- 
merical simulations show that event-chain algorithms clearly outperform the conventional Metropolis 
method. Irreversible versions of the algorithms, which violate detailed balance, improve the speed 
of the method even further. We also compare our method with a recent implementations of the 
molecular-dynamics algorithm. 



Hard spheres in three and in two dimensions (hard 
disks) occupy a special place in statistical mechanics. In- 
deed, many fundamental concepts, from the virial ex- 
pansion (by van der Waals and Boltzmann), to two- 
dimensional melting [l[, to long-time tails 0, were first 
discussed in these extraordinarily rich physical systems. 
These models have also played a crucial role in the his- 
tory of computation: both the Metropolis algorithm Q 
and molecular dynamics Q were first implemented for 
monodisperse hard disks in a box. In contrast with 
the spectacular algorithmic developments for lattice spin 
models Q , Monte Carlo algorithms for hard spheres 
have changed little since the 1950s, especially for high 
densities. Recent sophisticated implementation have re- 
duced the complexity of the molecular dynamics algo- 
rithm to a value comparable to that of the Monte Carlo 
method. Nevertheless, one can today still not equilibrate 
sufficiently large systems 0] to clarify whether the melt- 
ing transition in two-dimensional hard disks, at density 
(occupied volume fraction) 77 ~ 0.70, is weakly first or- 
der, or whether it is of the Kosterlitz-Thouless type [1], 
with a narrow hexatic phase in between the liquid and 
the solid. 

In this paper, we propose a class of Monte Carlo al- 
gorithms for hard-sphere systems: the "event-chain" al- 
gorithms. In contrast to the Metropolis algorithm, these 
methods are rejection- free. In a single move, they dis- 
place an arbitrary long chain of spheres, each advancing 
until it strikes the next one. Event-chain algorithms are 
generically faster than other Markov-chain algorithms, in 
part because the mean-square displacements of individ- 
ual particles are larger. In addition, one of the event- 
chain algorithms moves coherently over long distances. 
This further improves equilibration times. Finally, the 
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absence of rejections allows us to consider irreversible 
versions, which violate detailed balance, but preserve the 
correct stationary distribution. These versions accelerate 
the algorithm even further. The event-chain algorithms 
clearly outperform the traditional Metropolis algorithm 
for hard-disk and hard-sphere systems. 

In the local Metropolis algorithm, the move of a sphere 
is accepted if it induces no overlaps, and is rejected oth- 
erwise (see |Fig.T| ). This algorithm is notoriously slow at 
high density because, although a particle can move back 
and forth in the "cage" formed by its neighbors, it cannot 
easily escape from it 

To overcome the limitations of the local Metropolis 
algorithm, coordinated particle moves have been consid- 
ered: When the displacement of one sphere generates 
overlaps with other spheres, the latter are in turn moved 
out of the way. The rejection- free pivot cluster algorithm 
[loj . for example, works extremely well for binary 11] or 
for polydisperse [l^l mixtures, but it breaks down for 
the high densities of interest in two-dimensional melting. 
In Jaster's algorithm [l3| . overlapping spheres forming a 
chain are displaced, all of them by a fixed vector, until 
a configuration without overlaps is obtained (see |Fig.T| ). 
If a sphere branches out to more than one other sphere 
during the chain construction, the move is rejected (see 
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FIG. 1: Upper panels: Accepted (left) and rejected (right) 
local Metropolis moves of a disk in the cage formed by its 
neighbors. Lower panels: Accepted and rejected moves in 
Jaster's chain algorithm. 
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|Fig. This happens frequently, so the expected chain 
length is short and Jaster's algorithm barely faster than 
the local Metropolis algorithm. 




FIG. 2: Left two panels: Move of the straight event-chain 
(SEC) algorithm. The individual displacements add up to a 
distance £. Right two panels: Reflected event-chain (REC) 
move. 

In the algorithms presented here, each move consists in 
a deterministic chain of "events" : a disk advances until 
it strikes another one, which then in turn is displaced. 
The move starts with a randomly chosen disk, and stops 
when the lengths of all displacements add up to a total- 
displacement parameter £ (see |Fig.~2| ). This parameter 
allows the move to be reversible without rejections. To 
satisfy detailed balance, the move must also conserve 
configuration-space volume. This implies that when a 
disk strikes a neighbor, the latter may be displaced either 
in the original direction ("straight event-chain" (SEC) al- 
gorithm) or in the direction reflected with respect to the 
symmetry axis of the collision ("reflected event-chain" 
(REC) algorithm) (see |Fig.~2| ). In a periodic box, and 
with the initial direction 9 sampled uniformly in [0, 27r], 
both versions, which we call "reversible" , preserve the 
uniform measure because of detailed balance. 
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FIG. 3: Comparison of the integrated distribution of an ob- 
servable (the absolute value of the order parameter ^ of 
Eq. ([5])) between the SEC algorithm which breaks detailed 
balance, molecular dynamics (MD) and the local Metropolis 
algorithm for 1024 disks at 77 = 0.71. 



The detailed balance condition is allowed to be broken 
in the SEC algorithm. Indeed, for a given direction 9, a 
configuration F of disks can reach N other configura- 
tions in one move. By applying to F the N possible moves 
in direction —9, one finds the N configurations that reach 
F. Therefore, the SEC algorithm satisfies global bal- 



ance for any distribution of 9. Algorithms breaking de- 
tailed balance induce probability fiows in the configura- 
tion space and potentially speed-up equilibration [3l ■ We 
study such an irreversible version of the SEC algorithm 
where 9 is uniformly distributed in [0, tt]. To preserve er- 
godicity, at least two independent directions are needed. 
By far our fastest implementation (the "xy version" of 
the SEC algorithm) alternates moves in the positive x 
and y directions {9 = 0,7r/2). A version of the SEC 
algorithm, but with rejections and which cannot break 
detailed balance, was also mentioned in [isj . 

In |Fig. 3| we show the integrated distribution of l^*] of 
Eq. © 



^(l*l)d| *| 



(1) 



for the xy version of the SEC algorithm, for the Metropo- 
lis algorithm and for molecular dynamics in the same sys- 
tem. They are found to be equal. This demonstrates the 
correctness of our implementations. 

As a first step to analyze the performance of the 
event-chain algorithms, we consider the mean-square dis- 
placement (A^) of individual disks, both in the local 
Metropolis and in the event-chain algorithms. As men- 
tioned, event-chain algorithms generically outperform 
the Metropolis algorithm in part because they take larger 
steps on average. In order to compare the two methods, 
we measure time in units of attempted one-particle dis- 
placements. 




FIG. 4: Left: Histogram 7r(A/Ao) of the free path X(6) for 
1024 disks at density rj = 0.71. The distribution is close to 
exponential even in the solid phase. Right: Mean-free path 
Ao in units of the disk radius as a function of 77. 

Let us define the "free path" A = X{9) of a disk as 
the distance it must move in direction 9 to strike another 
particle. The ensemble average of A yields the mean-free 
path Aq. The distribution of the free path 7r(A/Ao) is well 
approximated by an exponential 



7r(A/Ao) ~ exp (-A/Aq) 



(2) 



even in the solid phase (see |Fig. 4[ ). This exponen- 
tial ansatz allows us to estimate the mean-square dis- 
placement for the local Metropolis algorithm, supposing. 
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for simplicity, that the proposed moves have fixed step 
size £ in random directions. The acceptance probability 
Pacc(£) = exp(-£/Ao) yields (A2(£)) = ee^^{-i/Xo), 



which is maximized when 



2Ao, 



:(A^(£)) = (A^(2Ao))=4A^/e2 



(3) 



To estimate the mean-square displacement for the 
event-chain algorithms, we suppose that the lengths of 
subsequent displacements in the chain are independent. 
In this case, the expected number of particles in the 
chain equals £/Ao -I- 1. We index the displacement dur- 
ing one event-chain move through a time-like variable s 
with < s < L The mean-square displacement of an 
event-chain move (the expected sum of the squares of 
the individual displacements) can be expressed through 
the probability 7r(s, s') that two times s and s' belong to 
the same particle movement: 

Jo Jo 

With the ansatz of Eq. dH), we have 7r(s,s') = 
exp (— |s — s'l/Ao). This yields the mean-square displace- 
ment per particle, which can be viewed as a short-time 
(local) diffusion coefficient: 



AocW 



This tends 



(A^(£)) 
{M{£)) 



= 2A; 



exp(-£/Ao) -I- i/Xo - 1 
i/Xo + 1 



(4) 



1 for l/Xo 



to 2Ag 



for large i/Xo, that is, to a value 
e^/2 ~ 4 times larger than what we found in Ec^. for 
the local Metropolis algorithm. This factor corresponds 
to the efficiency gain we might expect for a generic event- 
chain algorithm with large £/Xo, even though we will ob- 
tain considerably larger factors for the SEC algorithm. 

In a finite system, the expressions in Eq. (|4]) must be 
corrected for the center-of-mass displacement. For the 
SEC algorithm, the corrected mean-square displacement 
per particle, £'ioc(£), drops to zero for £/Ao oo because 
in that limit, for a finite box, all disks participate in the 
chain, and the system ends up moving as a solid block. 
The REC algorithm, in contrast, saturates to a constant 
mean-square displacement per particle for large chains. 

To further analyze the event-chain algorithms, we de- 
termined the auto-correlation time of the orientational 



order parameter \t [ij] for hard-disk systems at densities 
near the melting transition. The orientational order pa- 
rameter 4* averages the complex-valued local bond order 
ipj for each disk j, where 



and 



^ {k,j) 



(5) 



(6) 
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FIG. 5: Left: Order-parameter distribution for 256 disks in a 
periodic square box for rj = 0.71. Right: Correlation function 
C(At) for this system. The correlation time is obtained from 
an exponential fit, as shown. 



In Eq. ([S]), the sum is over the rij neighbors of j, and ipj^k 
is the angle of the shortest vector equivalent to — Xj 
14j . Probable values of ^ form an irregular ring around 
the origin (see the scatter plot in |Fig. 5[ the 



* + TT 

symmetry in a square box imposes (^'(t)) ~ 0). 

We suppose that the correlation time of this system 
is proportional to the time the order parameter takes 
to wander around the ring, that is, the auto-correlation 
time of ^. This global measure of the overall speed of 
a Monte Carlo simulation is more appropriate than, for 
example, single-particle diffusion constants. However 
is very long to decorrelate at the interesting densities (see 
|Fig. 7[ ), and we have to limit ourselves to small systems. 
The auto-correlation function C{At) of the orientational 
order parameter is defined by the ensemble average 

normalized so that C(0) = 1. In our square box, this 
function decays to zero exponentially for large At (see 
|Fig. 5[ ), and the decay constant r and the speed l/r 
arc obtained by a fit, for each value of the parame- 
ters {N,7^,t}, from one single very long simulation (with 
t ^ r). The local Metropolis algorithm, for its optimal 
step size, is as fast as the event-chain algorithms with 
1/Xq ~ 1. (Our implementation moves 3 x 10^° particles 
per hour on a 2.8GHz single-processor workstation for 
N = 4096.) 

For small total displacements £/Ao <C 1, the speed 
of all the algorithms (reversible and irreversible SEC, 
REC, and local Metropolis algorithm) is proportional 
to D^^{tj, as given by Eq. that is, they aU follow 
the single-particle behavior (see |Fig."6l ). For larger l/Xo, 
the event-chain algorithms realize a considerable speed- 
up compared to the local Metropolis algorithm (also in 
|Fig. 6[ ). Moreover, both versions of the SEC algorithm set 
up coherent motion across the system and are clearly bet- 
ter than the REC algorithm, whose particle chains mean- 
der through the system (as shown in |Fig. 2| , so that the 
disks move incoherently. For large ^/Aq, the irreversible 
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FIG. 6: Left: Efficiency of SEC and REC algorithms for 1024 
disks at ?7 = 0.71 (all speeds normalized by the speed of the 
reversible SEC algorithm at ^/Ao = 1). The speed of the 
local Metropolis algorithm and the mean-square displacement 
per particle from Eq. Q are also shown. Right: Density 
dependence of the optimal speed-up factor. 



SEC algorithm is faster than the reversible version: it is 
of advantage to break detailed balance. [Figure 6| also il- 
lustrates that the SEC algorithm becomes more efficient 
(as compared to the local Metropolis algorithm) as one 
approaches the transition from the liquid phase (at den- 
sity 77 ~ 0.708). The optimal speed-up increases with the 
system size, as shown in lTableTl This suggests that the 
speed-up of the SEC algorithm may well increase with 
the correlation length of the system, and may, in the 
transition region, have a more favorable scaling than the 
local Metropolis algorithm. 





Optimal speed-up 


iV 


Reversible 


Irreversible 


64 




8 


256 


~ 8 


~ 11 


1024 


~ 9 


~ 15 


4096 


~ 10 


~ 20 



TABLE I: Optimal speed-up reached by the SEC algorithm 
(with respect to the reversible SEC algorithm for £/Xo = 1) 
at density rj = 0.71 as a function of particle number. 

Let us finally discuss the relationship between the 
Monte Carlo method and the molecular-dynamics algo- 
rithm. All these approaches describe the same equilib- 
rium state. Unlike the Monte Carlo method, the molec- 
ular dynamics follows the physical time-evolution of the 
system. The first implementations of the molecular dy- 
namics algorithm Q were very time consuming, with a 
complexity of 0{N) per event (collision), slower than the 
Metropolis algorithm (0(1) per move). The complexity 
of modern implementations has improved to O(logA^) 
15 1 per event and even 0(1) [l^]. A closer look is thus 
needed to choose between the two methods. 

We used a simple version of the molecular dynamics to 



compute the decorrelation time of in the same way as in 
|Fig. 5| In number of events, molecular dynamics is found 
to be about three times faster than the irreversible ver- 
sion of SEC for 77 ~ 0.7 and iV = 64 - 1024. It is very in- 
teresting to notice that molecular dynamics shows, unlike 
REC, the same density dependence of its speed as SEC 
around the transition region. We then determined the 
CPU time per collision of one of the most rapid current 
implementations of the hard-disk molecular-dynamics al- 
gorithm For the 32 x 32 system at 77 ^ 0.7, this 
implementation reaches about 1.7 x 10^ collisions per 
hour on a 2.6GHz workstation [l7j . Our xy implemen- 
tation of the SEC algorithm reaches 3 x 10^° collisions 
per hour on similar hardware. Our implementation is 
thus about 5 times faster in CPU-time to reach thermo- 
dynamic equilibrium than the best molecular-dynamics 
implementation. We should also note that SEC is much 
easier to implement. A synopsis of these relative and ab- 
solute timing issues is presented in |Fig. 7| For clarity, 
we give times in terms of "equivalent Metropolis moves" , 
this means that one event of the molecular dynamics al- 
gorithm corresponds to '--^ 3 SEC events and to ^ 60 
Metropolis moves. The horizontal lines indicate what can 
be achieved in approximately one hour with our imple- 
mentation of the Metropolis algorithm, irreversible SEC, 
and the implementation of the molecular dynamics algo- 
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FIG. 7: System-size dependence of the correlation time of the 
orientational order parameter for two densities. What can 
be achieved in approximately one hour using the different 
algorithms discussed in this paper is indicated by horizontal 
lines. 

In conclusion, we have in this paper proposed a class of 
algorithms for hard spheres and related systems, which 
clearly outperform the local Metropolis algorithm. We 
discussed three aspects of our algorithms, which all con- 
tributed to improve their speed. First, we showed that 
event-chain algorithms have a larger effective step size 
than the local Metropolis algorithm, because spheres 
move until they strike one of their neighbors. We com- 
puted mean-square displacements per particle (local dif- 



5 



fusion constants) to quantify this point. Nevertheless, 
local diffusion constants are not clearly related to the 
speed of the algorithm: they merely describe the short- 
time rattling of a particle in its cage (only for small ^/Aq 
is the local diffusion constant directly proportional to the 
algorithm's speed). Second, we performed numerical sim- 
ulations of two variants of the method, and carefully an- 
alyzed the auto-correlation function of the orientational 
order parameter. One of them, the SEC algorithm, in- 
duces coherent motion of a long chain of spheres, and it 
allows the different parts of the system to communicate 
with each other. We witnessed considerable performance 
gains of this algorithm in the critical region, in the same 
way than the molecular dynamics. This suggests the ex- 
citing possibility that the speed-up of the event-chain al- 
gorithm grows with the correlation length of the system, 
and may have a more favorable scaling than the local 
Metropolis algorithm in the critical region. This speed- 
up, which is shared by both the molecular-dynamics al- 
gorithm and the SEC algorithm, is not understood and 
should be further investigated. Third, we noticed that 
the absence of rejections permitted to conceive an irre- 
versible version of the SEC algorithm which improves the 
performances. 

Our implementation of the SEC algorithm approaches 
equilibrium (for large systems at 77 ~ 0.70) about 40 
times faster than our local Metropolis algorithm, not 
only because of the speed-up evidenced in |Fig. 6| but 
also because the xy version of the algorithm computes 
no scalar products and uses very few random numbers. 
It also equilibrates about five times faster than the best 
molecular-dynamics implementation and preserves cer- 
tainly a large potential for improvement. 

Nevertheless, CPU times needed for convergence re- 
main extremely large, and even with our algorithm, full 
convergence of systems with 10^ particles at high densi- 
ties comes barely into reach. The irreversible SEC algo- 
rithm not only appears to be the fastest currently known 
simulation method for dense hard-disk and hard-sphere 
systems, but it also provides a telling example of the 
benefits of breaking detailed balance in Monte Carlo al- 
gorithms going beyond the "lifting" Markov chains [l^ . 
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